2 Radar Data Preprocessing

Weather radar data of the firework events at the turns of the years usually contain some degree of precipitation clutter. To filter out precipitation advanced algorithms such as MistNet have been developed, but as we are dealing with dual-pol data here, we can use a simpler and yet robust method using the depolarization ratio (Kilambi, Fabry, and Meunier 2018).

To make sure our processed weather radar data does not contain any significant proportions of precipitation or ground clutter anymore, we process the data as follows:

  1. We remove electromagnetic interference based on a visual inspection of the scan and throw out all data of affected rays.
  2. We calculate the depolarization ratio (Kilambi, Fabry, and Meunier 2018) and separate biology from meteorology by classifying all range gates with a depolarization ratio \(>-12dB\) as biology. We subsequently ‘despeckle’ this, to remove obvious misclassifications.
  3. We average reflectivity over a number of scans before the time of the fireworks event and throw out the range-gates with highest average reflectivities.

All these steps can be undertaken directly on the polar volume data, so we can subsequently directly plug the cleaned up volume into the range-bias correction.

2.1 Setting-up the pre-processing environment

As usual, we use bioRad (Dokter et al. 2019), but this time we include plotly for some interactive plotting.

2.2 Removing electromagnetic interference

We have determined in which scans birds are taking off based on the maximum increase in reflectivity in the scan for each of the involved radars. Let’s now look at these scans to see how much filtering for electromagnetic interference we need to do. The easiest way to determine which rays are subject to this interference is by plotting the scans in polar coordinates (\((r, \alpha)\)), so interference stands out as horizontal lines of more or less constant, or very gradually changing reflectivities. Plotting using plotly makes it easier to identify the specific problematic rays as one can zoom in to identify the exact azimuths \(\alpha\).

The scans we will be using:

  • Herwijnen: RAD_NL62_VOL_NA_201712312310_ODIM.h5
  • Den Helder: RAD_NL61_VOL_NA_201712312305_ODIM.h5

For illustrative purposes we will only illustrate removal of EM interference for the Herwijnen, as the procedure for Den Helder is exactly identical, but this scan contains very little of said clutter.

Right away we can see that rays at two places in the scan are subject to electromagnetic interference. This is probably most problematic in the lowest elevations of the volume scan, but nevertheless each of the length(pvol_hrw$scans) scans have be checked manually. Doing so results in the identification of the following rays that contain electromagnetic interference (ei_rays), organised in a list with the scan numbers (organised ascendingly per elevation angle) as keys. Admittedly: there is another ray that seems to contain interference in the first scan, but this is so far away from the radar (150km+) it should not affect our results as no meaningful numbers of birds can be detected at that range anyways. Similarly, there are similar patterns of interference/clutter in higher elevation scans, but these too should not affect our results.

We can now remove the data for the affected rays in the corresponding scans by setting the values to NA (see R/remove_rays.R).

2.2.1 Verify removal of rays with EM interference

If removal is correct, the \((r,\alpha)\) plots should not show clear horizontal structures anymore.

Well, that seems to work nicely. The remaining bits of clutter won’t affect the results much as they cover fairly small areas and consistently appear at longer distances away from the radar, where birds are difficult to detect anyways.

2.3 Filter meteorology using the depolarization ratio

Meteorology can be filtered using the depolarization ratio following Kilambi et al. (2018). We calculate the depolarization ratio for the raw pvol data after EM interference has been removed and subsequently ‘despeckle’ the results to improve the classification.

Despeckling works by comparing the classification of the majority of the neighbourhood rangegates with the classification of the center rangegate, and changing the latter to reflect the majority of the neighbourhood classification if there is a difference. We define the ‘neighbourhood’ as a \(3^{\circ}\) by \(3 \times rscale\) area centered around a focal rangegate (3 rangegates in azimuth \(\times\) 3 rangegates in range) . Selecting the rangegates while taking the sphericity of the radar scan into account (e.g. ray 360 should be directly adjacent to ray 1) is made easier with the R/window_coords.R function. The despeckling is implemented in R/despeckle_scan_logical.R.

With the despeckling algorithm in place, we can:

  1. Calculate the depolarization ratio (DPR).
  2. Classify biology as rangegates where DPR > -12 and store this classification as BIOLR (Biology Raw) scan parameter in the pvol object.
  3. Despeckle the classification and store the outcome in the BIOLD (Biology Despeckled) scan parameter in the pvol object.
source("R/window_coords.R")
source("R/despeckle_scan_logical.R")

# Calculate depolarization ratio, classify and despeckle biology classifications for the entire volume
calculate_dpr <- function(pvol){
  for (i in seq_along(pvol$scans)) {
    # Calculate ZDR as ZDR = DBZH - DBZV
    pvol$scans[[i]]$params$ZDR <- pvol$scans[[i]]$params$DBZH - pvol$scans[[i]]$params$DBZV
    attributes(pvol$scans[[i]]$params$ZDR)$param <- "ZDR"
    
    # Calculate depolarization ratio
    zdr_linear <- 10 ** (pvol$scans[[i]]$params$ZDR / 10)
    dpr_linear <- (zdr_linear + 1 - 2 * sqrt(zdr_linear) * pvol$scans[[i]]$params$RHOHV) / 
                  (zdr_linear + 1 + 2 * sqrt(zdr_linear) * pvol$scans[[i]]$params$RHOHV)
    pvol$scans[[i]]$params$DPR <- 10 * log10(dpr_linear)
    attributes(pvol$scans[[i]]$params$DPR)$param <- "DPR"
    
    # Classify based on depolarization ratio
    biology <- (pvol$scans[[i]]$params$DPR > -12) * 1  # multiply by 1 to convert TRUE/FALSE to 1/0
    class(biology) <- c("param", "matrix")
    attributes(biology) <- attributes(pvol$scans[[i]]$params$DPR)  # copy attributes from DPR
    attributes(biology)$param <- "BIOLR"
    pvol$scans[[i]]$params$BIOLR <- biology
    
    # Despeckle biology classification
    pvol$scans[[i]]$params$BIOLD <- pvol$scans[[i]]$params$BIOLR
    pvol$scans[[i]]$params$BIOLD <- despeckle_scan_logical(pvol$scans[[i]]$params$BIOLD)
    attributes(pvol$scans[[i]]$params$BIOLD)$param <- "BIOLD"
    
    # Classify based on correlation coefficient
    biology_rhohv <- (pvol$scans[[i]]$params$RHOHV < 0.95) * 1
    class(biology_rhohv) <- c("param", "matrix")
    attributes(biology_rhohv) <- attributes(pvol$scans[[i]]$params$RHOHV)
    attributes(biology_rhohv)$param <- "BIOLRRHOHV"
    pvol$scans[[i]]$params$BIOLRRHOHV <- biology_rhohv
    
    # Despeckle biology classification based on correlation coefficient 
    pvol$scans[[i]]$params$BIOLDRHOHV <- pvol$scans[[i]]$params$BIOLRRHOHV
    pvol$scans[[i]]$params$BIOLDRHOHV <- despeckle_scan_logical(pvol$scans[[i]]$params$BIOLDRHOHV)
    attributes(pvol$scans[[i]]$params$BIOLDRHOHV)$param <- "BIOLDRHOHV"
  }
  
  return(pvol)
}

pvol_hrw <- suppressWarnings(calculate_dpr(pvol_hrw)) # Will throw NaN warnings if not suppressed
pvol_dhl <- suppressWarnings(calculate_dpr(pvol_dhl))

2.3.1 Verify DPR-based classification

Now let’s plot some PPIs to verify the accuracy of DPR-based classification and the subsequent despeckling, by plotting DBZH, VRADH, DPR, BIOLR and BIOLD.

The plots show accurate classification of the obvious precipitation zones, except at the edges of these echoes, where BIOLD is a vast improvement over BIOLR, showing the value of despeckling. Similarly, there is a lot of ‘noise’ where birds should be, but despeckling takes care of most of that quite nicely as well. Additionally, it shows a pattern we would expect to see: at closer distances to the radar most ‘speckles’ that are not near to precipitation zones are turned into biology, and at distances further from the radar they are more often ‘flipped’ to meteorology. This method may not be perfect, but it classifies birds quite conservatively. The few misclassifications that remain should not affect the results so much, as they are few in number and do not occur at the centers of precipitation echoes, so they are not likely to turn into numerical outliers.

2.4 Remove classified precipitation from polar volumes

Now that we have accurate classifications of the rangegates based on depolarization ratios, we can start to remove the precipitation from the polar volumes, to retain a scan that comprises of only birds (with a few occasional misclassifications). As there are areas where DPR and DBZH do not overlap, we also have to remove all rangegates that are not classified.

Plotting the same PPIs as before should now show a cleaned-up/precipitation-free scan next to the classifications.

That looks very good for both Herwijnen and Den Helder radars, but for the latter we have a lot of sea clutter that still needs to be removed, but that will come when ground clutter is filtered.

2.5 Filter ground clutter

We will filter out ground clutter by calculating summary statistics of the rangegate reflectivities over:

  1. The 36 scans preceding the scans selected for the study of the fireworks event (= 3 hours worth of scans).
  2. A day of clear weather closest to the 31st of December 2017.

For each we will filter ground clutter based on the mean DBZH values. Using the variance and mad of the DBZH was tested, but has a few difficulties:

  1. variance is very sensitive to the outliers caused by rangegates with NA values (detection below the ‘mds’, the minimum detectable signal) occasionally flipping over to a noisy measurement, resulting in very high variances.
  2. mad is much more robust to outliers, but to compute these values we need to set NA cells to the ‘mds’ (minimum detectable signal), which will result in mad values close to, or exactly 0 for cells that never reflected as well as true static clutter, so it’s difficult to separate those.

Finally, a visual inspection showed the mean and mad of DBZH (assuming one could overcome the aforementioned problem with the latter) do not differ much, but the mean is somewhat more ‘aggressive’ in filtering, which in this case is quite good.

Combining the clutter removal based on a clear day as well as the 36 preceding scans lets us account for both truly static clutter (e.g. buildings) as well as clutter that is more dynamic such as sea and wind park clutter, without also requiring us to resort to filtering of dynamic clutter using a VRADH threshold. The quality of filtering is assessed visually.

2.5.1 Dynamic clutter

We select 36 (3 hours worth of scans) preceding the start of the fireworks (23:00 UTC) and add an additional margin of 3 scans (15 minutes of scans) as that the VIR plots in the previous chapter have shown numbers of birds aloft are very low and stable up to that period.

We can now loop over the files one by one and stack reflectivity data (DBZH) — after filtering out precipitation — in a multidimensional array.

Note: the following code chunk will only run in full-reproduction mode as it takes quite a lot of time. Results are saved, so the next iteration this chunk can be skipped.

With all DBZH compiled in a single multidimensional array, we can calculate mean reflectivity, which we store as DBZH_AVG in a pvol that now contains the dynamic clutter map.

2.5.1.1 Verify dynamic clutter map

Let’s see what that looks like on a basemap, using a DBZH_AVG threshold of \(-10dbZ\), following (Dokter et al. 2011).

Visually assessing this clutter map shows that it works quite well, selecting e.g. areas with wind parks, sea clutter, high buildings, industry, etc. Exactly what we hoped to achieve.

2.5.2 Static clutter

Now, let’s retry exactly the same procedure, but this time selecting a day with no precipitation, which can be done using this tool by KNMI, so we can filter for truly static clutter.

We select the following days:

  1. Herwijnen: December 29th, 2017
  2. Den Helder: December 25th, 2017

Note: the following code chunk will only run in full-reproduction mode as it takes a lot of time to run. Results are saved, so the next iteration this chunk can be skipped.

And we calculate mean DBZH values (DBZH_AVG).

2.6 Range-bias correction

With all identifiable sources of clutter removed from the raw polar volume, we can apply the range-bias correction (Kranstauber et al. 2020). For this it is necessary to calculate the local vertical profile for each of the radars. Ideally, this would be done using the filtered pvol we have now generated, but the vol2bird algorithm (Dokter et al. 2011) only takes pvol files as input, rather than R objects. As there is no implementation of a converter yet, for now a vp of the raw pvol files will have to do. As there is no precipitation within the relevant distance to the radars (5-35km), the calculated vp based on the raw pvol files should not differ wildly from that of the filtered pvol R object we have generated in the previous steps.

For the Den Helder radar we calculate the vp by setting azimuthal limits to cover the mainland of North Holland, rather than the whole radar domain, as the latter will result in vps that underestimate the true density of birds aloft. See Appendix [generating-vps-for-den-helder-radar] for a more detailed explanation.

We can now plot the final range-corrected PPIs.

It’s a little hard to interpret with the PPI pixels that contain no birds set to 0 and no landscape below, so let’s see what it looks like on a basemap. As there are some ships out at the North Sea causing reflectivities orders of magnitude higher and thus stretching the colormap, we — for now — need to ‘clip’ these values to a maximum set at the value of the 99th percentile.

2.6.1 Range-corrected PPI of Herwijnen radar

2.6.2 Range-corrected PPI of Den Helder radar

References

Dokter, Adriaan M., Peter Desmet, Jurriaan H. Spaaks, Stijn van Hoey, Lourens Veen, Liesbeth Verlinden, Cecilia Nilsson, et al. 2019. “bioRad: Biological Analysis and Visualization of Weather Radar Data.” Ecography 42 (5): 852–60. https://doi.org/10.1111/ecog.04028.

Dokter, Adriaan M., Felix Liechti, Herbert Stark, Laurent Delobbe, Pierre Tabary, and Iwan Holleman. 2011. “Bird Migration Flight Altitudes Studied by a Network of Operational Weather Radars.” Journal of the Royal Society Interface 8 (54): 30–43. https://doi.org/10.1098/rsif.2010.0116.

Kilambi, Alamelu, Frédéric Fabry, and Véronique Meunier. 2018. “A Simple and Effective Method for Separating Meteorological from Nonmeteorological Targets Using Dual-Polarization Data.” Journal of Atmospheric and Oceanic Technology 35 (7): 1415–24. https://doi.org/10.1175/JTECH-D-17-0175.1.

Kranstauber, Bart, Willem Bouten, Hidde Leijnse, Berend-Christiaan Wijers, Liesbeth Verlinden, Judy Shamoun-Baranes, and Adriaan M Dokter. 2020. “High-Resolution Spatial Distribution of Bird Movements Estimated from a Weather Radar Network.” Remote Sensing 12 (4). Multidisciplinary Digital Publishing Institute: 635.